rm(list = ls())
cat("\014")

library(tidyverse, quietly = T)
## Warning: package 'tidyverse' was built under R version 3.5.1
## -- Attaching packages ---------------------------------------------------------------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 3.1.0     v purrr   0.2.5
## v tibble  2.0.1     v dplyr   0.7.8
## v tidyr   0.8.1     v stringr 1.3.1
## v readr   1.3.1     v forcats 0.3.0
## Warning: package 'ggplot2' was built under R version 3.5.2
## Warning: package 'tibble' was built under R version 3.5.2
## Warning: package 'tidyr' was built under R version 3.5.1
## Warning: package 'readr' was built under R version 3.5.2
## Warning: package 'purrr' was built under R version 3.5.1
## Warning: package 'dplyr' was built under R version 3.5.2
## Warning: package 'stringr' was built under R version 3.5.1
## Warning: package 'forcats' was built under R version 3.5.1
## -- Conflicts ------------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(sf,        quietly = T)
## Warning: package 'sf' was built under R version 3.5.2
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(gganimate, quietly = T)
library(openair,   quietly = T)
## Warning: package 'openair' was built under R version 3.5.1
library(lubridate, quietly = T)
## Warning: package 'lubridate' was built under R version 3.5.1
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
library(zoo,       quietly = T)
## Warning: package 'zoo' was built under R version 3.5.1
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(ggthemes,  quietly = T)
## Warning: package 'ggthemes' was built under R version 3.5.2
library(ggmap,     quietly = T)
## Warning: package 'ggmap' was built under R version 3.5.2
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library(sp,        quietly = T)
## Warning: package 'sp' was built under R version 3.5.1
library(magick,    quietly = T)
## Warning: package 'magick' was built under R version 3.5.2
## Linking to ImageMagick 6.9.9.14
## Enabled features: cairo, freetype, fftw, ghostscript, lcms, pango, rsvg, webp
## Disabled features: fontconfig, x11

Import the time-activity diary for the line and tidy it up.

diary      <- read_csv('air_quality_data/victoria_line_2014-12-08.csv', col_types = cols())

base       <- tibble(date_time = seq(min(diary$start_time), max(diary$end_time), by = 'min'))

base       <- left_join(base, diary, by = c('date_time' = 'start_time')) %>% 
              select(-fake_id)

base       <- left_join(select(base, -end_time),
                        filter(base, date_time != end_time),
                        by = c('date_time' = 'end_time')) %>%
              mutate(line.x            = if_else(is.na(line.x)            & !is.na(line.y),            line.y,            line.x),
                     station.x         = if_else(is.na(station.x)         & !is.na(station.y),         station.y,         station.x),
                     environment.x     = if_else(is.na(environment.x)     & !is.na(environment.y),     environment.y,     environment.x),
                     sub_environment.x = if_else(is.na(sub_environment.x) & !is.na(sub_environment.y), sub_environment.y, sub_environment.x)) %>%
              select(date_time, line.x, station.x, environment.x, sub_environment.x) %>%
              rename(line = line.x, station = station.x, environment = environment.x, sub_environment = sub_environment.x)
## Warning: package 'bindrcpp' was built under R version 3.5.1
rm(diary)

Import the measured air quality to join to the time-activity diary

air_quality <- read_csv('air_quality_data/LUexportMay2016.csv', col_types = cols()) %>%
                rename_all(tolower) %>%
                mutate(datetime = as.POSIXct(datetime, format='%d/%m/%Y %H:%M')) %>%
                filter(sitecode == 'CAR')

Join the time activity and air quality

line       <-  left_join(base, air_quality, by = c('date_time' = 'datetime')) %>% 
                mutate(scaledvalue = if_else(species == 'PM25', value, scaledvalue)) %>%
                select(-value) %>%
                rename(concentration = scaledvalue)

rm(air_quality, base)

Import the background air quality for that day to correct the PM2.5 measurements

background_pm25 <- importKCL(site = "kc1", year = c(2014,2015,2016), pollutant = "pm25", met = FALSE,
                             units = "mass", extra = FALSE)
## NOTE - mass units are used 
## ug/m3 for NOx, NO2, SO2, O3; mg/m3 for CO
## PM10_raw is raw data multiplied by 1.3
background_pm25 <- data.frame(background_pm25, day = as.Date(format(background_pm25$date)))
background_pm25 <- aggregate(pm25 ~ day, background_pm25, mean)

Join the background PM2.5 to the monitored concentrations

line     <-    mutate(line, day = date(date_time)) %>%
                left_join(background_pm25, by = c('day' = 'day')) %>%
                select(-day)

Correct the PM2.5 data

line <-         mutate(line, concentration = 
                         case_when(concentration >  pm25/0.44 & species == 'PM25' & !is.na(concentration) ~ pm25 + (1.82 * (concentration - (pm25/0.44))),
                                   concentration <= pm25/0.44 & species == 'PM25' & !is.na(concentration) ~ concentration * 0.44,
                                   TRUE ~ concentration)) %>%
                select(-pm25)

Get the station locations

stations <- st_read('https://raw.githubusercontent.com/dracos/underground-live-map/master/bin/stations.kml', quiet = T) %>%
            select(-Description) %>%
            rename_all(tolower) %>%
            mutate(name = gsub(' Station', '', name)) %>%
            mutate(name = gsub("'", '', name)) %>%
            st_zm(drop = TRUE)

Join locations to the concentrations and station names

line   <-left_join(line, stations, by = c('station' = 'name')) %>% st_as_sf()

Fill in the missing coordinates using linear interpolation

st_geometry(line) <- as_tibble(st_coordinates(line)) %>% 
                      mutate(X = na.approx(X), Y = na.approx(Y)) %>% 
                      st_as_sf(coords = c('X', 'Y')) %>%
                      st_geometry()

Start trying to animate. First the graph animation.

text_placement <- filter(line,species == 'PM25') %>%
                    as_tibble() %>%
                    select(-geometry) %>%
                    group_by(species) %>%
                    summarise(max_x = max(date_time), max_y = max(concentration))

line           <- left_join(line, text_placement, by = c('species' = 'species'))

graph <- ggplot(data = filter(line,species == 'PM25')) +
    geom_path(aes(date_time, concentration), colour = '#0099CC', size = 1.2) +
    geom_point(aes(date_time, concentration), colour = 'red', size = 4) +
    geom_text(aes(max_x, max_y,label = round(concentration,0)),
              size = 10, hjust = 1, vjust = 1, colour = '#0099CC') + 
    theme(axis.text.x  = element_blank(),
          axis.text.y  = element_text(size = 12),
          axis.title.y = element_text(size = 12),
          axis.line    = element_line(colour = 'black')) +
    xlab('') +
    ylab('PM2.5 ug/m3') +
    transition_reveal(date_time)

graph_animation <- animate(graph, fps = 3, height = 600)
anim_save(file='graph_animation.gif', animation = graph_animation, path = 'outputs')

Now the spatial animation.

line <- st_set_crs(line,4326)

map_area <- c(as.vector(st_bbox(line))[1]-0.01, as.vector(st_bbox(line))[2]-0.01, as.vector(st_bbox(line))[3]+0.01, as.vector(st_bbox(line))[4]+0.01)

map <-  ggmap(get_map(location = map_area), source = "google", maptype = "roadmap", zoom = 13) +
    geom_path(data = as_tibble(st_coordinates(filter(line,species == 'PM25')[1:56,])), aes(X, Y), colour = "#0099CC", size = 1) +
     geom_sf(data = filter(line,species == 'PM25'), size=4, colour = 'red',inherit.aes = FALSE) +
     coord_sf(datum=NA)  +
     transition_time(date_time) +
     ease_aes('linear') +
     theme(axis.title = element_blank(),
           panel.background = element_blank(),
           panel.border = element_rect(colour = 'black', fill = NA),
           axis.text = element_blank())
## Source : http://tile.stamen.com/terrain/13/4092/2720.png
## Source : http://tile.stamen.com/terrain/13/4093/2720.png
## Source : http://tile.stamen.com/terrain/13/4094/2720.png
## Source : http://tile.stamen.com/terrain/13/4092/2721.png
## Source : http://tile.stamen.com/terrain/13/4093/2721.png
## Source : http://tile.stamen.com/terrain/13/4094/2721.png
## Source : http://tile.stamen.com/terrain/13/4092/2722.png
## Source : http://tile.stamen.com/terrain/13/4093/2722.png
## Source : http://tile.stamen.com/terrain/13/4094/2722.png
## Source : http://tile.stamen.com/terrain/13/4092/2723.png
## Source : http://tile.stamen.com/terrain/13/4093/2723.png
## Source : http://tile.stamen.com/terrain/13/4094/2723.png
## Source : http://tile.stamen.com/terrain/13/4092/2724.png
## Source : http://tile.stamen.com/terrain/13/4093/2724.png
## Source : http://tile.stamen.com/terrain/13/4094/2724.png
## Source : http://tile.stamen.com/terrain/13/4092/2725.png
## Source : http://tile.stamen.com/terrain/13/4093/2725.png
## Source : http://tile.stamen.com/terrain/13/4094/2725.png
## Source : http://tile.stamen.com/terrain/13/4092/2726.png
## Source : http://tile.stamen.com/terrain/13/4093/2726.png
## Source : http://tile.stamen.com/terrain/13/4094/2726.png
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
map_animation  <- animate(map, fps = 3, height = 600)

anim_save(file='map_animation.gif', animation = map_animation, path = 'outputs')

Print

graph_animation

map_animation